ROCKi
1 Samples
2 Check for KO of the MIC
MIC inclusion:
Note
Two points from the WTs, one in each group, are not plotted since the junction coverage is not enough to calculate PSI.
3 PCA
In order to consider a gene for further analysis, it must be expressed >20 reads in all 3 replicates from the at least a group. Points are colored based on the sample group, as the PSI plot.
As it can be appreciated it is not clear the groups cluster together.
“I guess, we can also say that PC2 loosely and dirtily separates ROCKi vs Veh” (Manu, 15/08/2023)
I would rather say that PC2 kind of separates KO Veh from the rest of the samples, but it is not that clear:
4 Multiple comparisons
As in a control-genotype + treatment experiment many comparisons can be calculated. Let’s first start with WT + Veh as the reference and identify how many genes are differentially expressed. For comparisons of two groups IHW was used for correcting p-values and apeglm for log2fold change shrinkage. Genes were considered differentially expressed if corrected p-value <0.1 disregarding the corrected log2foldchange.
4.0.1 KO Veh vs WT Veh
Only 37 genes differentially expressed
4.0.2 WT ROCKi vs WT Veh
As identified by Manu, there are no big differences, here I just could identify a single gene:
4.0.3 KO ROCKi vs KO Veh
Here the reference changes and now is the KO treated with Veh, this was done to identify genes that their transcription rate changed with the ROCKi. Unfortunatly only one showed statistical differences.
5 KO changes
We could then use the samples to identify genes which expression is change in the KO disregarding the treatment. For this we used a model ~ Genotype + Treatment and looked for differences between WT and KO. Using same parameters as before I get 134 genes, which are plotted in the heatmap.
5.1 GO enrichment
5.1.1 Biological Process
Bubbles represent significant enriched GO terms, alpha and size show significance. Terms are grouped with Wang distance into bubbles unfilled. Most significant term is shown as label.
5.1.2 Molecular Function
5.2 KEGG enrichment
6 Treatment changes
Similar as with the KO changes I tested then if disrigarding the genotype we could identify transcriptional changes related to the ROCKi treatment. As expected from the comparison WT ROCKi vs WT Veh not many genes showed differential expression, only one.
In order to explore further ROCKi, here is the KEGG pathway “Regulation of actin cytoskeleton” where ROCK plays an important role. It is wroth to mention that ROCK acts at the level of proteins,phosphorylating its targets. As appreciated no changes are observed!!!
7 Merging comparisons.
With Patryk’s mail I got inspired to look for specific patterns,as the one he was pointing out, which are not being recovered in the previous analysis. I thought with the differential expression analysis I could get there, but with the stringent (I tried more relaxed ones with similar results) parameters I used, very very very few genes were found differentially expressed and almost no overlap between lists (if any). Looking at the raw p-values I realized that some genes presented low p-values for multiple comparisons but when adjusting them significance is lost. Since I was looking for a specific pattern where multiple comparisons need to be combined I found a couple of ways for summarizing multiple p-values (from independent tests) into one (more info here).
Patterns to look for (NOT FILTERS):
KO mimicking; ROCKi shifts WT expression in Veh into KO.
- WT vs KO should be different treated with Veh
- Veh vs ROCKi should be equal in KO ( therefore WT Veh vs KO ROCKi should be different) **
- Veh vs ROCKi should be different in WT (same direction as WT vs KO)
KO rescue; ROCKi shifts KO expression in Veh into WT.
- WT vs KO should be different treated with Veh
- Veh vs ROCKi should be different in KO (same direction as WT vs KO)
- Veh vs ROCKi should be equal in WT ( therefore KO Veh vs WT ROCKi should be different) **
Important
** Since we are looking for small p-values I will use the comparisons where differential expression is looked for. The comparison where “should be equal” will not be considered and p-value ignored.
Genes which expression is altered by KO disregarding treatment were not considered in this part, because result would be then biased through genotype differences.
Genes were considered for this part if in all comparisons at least all samples from a single group showed experession (<20 counts).
Note
Not all comparisons have the same importance, therefore Stouffer method was used. With this method it is common to assign weights to p-values. In this part I used a ratio 1:1:2 for weights having bigger weights the comparison Veh vs ROCKi in WT, for KO mimicking, and in KO, for KO rescue. The result of this method gives a new p-value for each gene which was then corrected with IHW, using as covariates the mean of TPM.
7.1 KO mimicking
Comparisons Names Codes:
kovswt-Veh –> KO Veh vs WT Veh
rvsv-WT –> WT ROCKi vs WT Veh
rvsv-KOWT –> KO ROCKi vs WT Veh
Only one gene showed this pattern of expression.
Caution
Patryk’s example Nr4a2 showed this pattern, unfortunately I did not found it with this approach. Nevertheless, when I looked at the expression of that same gene I found that I did not get the same gene counts. I don’t know why is that so different since the others rest look quite similar. (I did not considered expressed Abcc6 and Lrmp ->Irag2)
7.2 KO rescue
Comparisons Names Codes:
kovswt-Veh –> KO Veh vs WT Veh
rvsv-KO –> KO ROCKi vs KO Veh
rvsv-WTKO –> WT ROCKi vs KO Veh
Important
kovswt-Veh log2FoldChanges should be interpreted here as the opposite since we are focusing on the changes realted to KO Veh.
7.2.1 GO enrichment (BP)
Here it is important to mention that the universe is reduced since not all expressed genes were tested. White boxes represent enriched go terms, size and alpha is scaled based on significance. Additionally, GO terms are grouped in black boxes, based on “wang” distances. Most significant term from cluster is shown as parent “term”.
7.2.2 GO enrichment (MF)
7.2.3 KEGG enrichment
7.3 LTR
THe likelihood ratio test from DESeq2 is an alternative to pair-wise comparisons to analyze all levels of all factors at once. When implemented, with IWH for p-value correction, we could identify 62 differentially expressed genes. Which expression looks like this:
55 of these belong to genes previously characterized as KO changes which pattern is recognizable in the heatmap. The rest look like this:
Two of them have a pattern of KO rescue (Gm19680 and Ntng2). Fam205a2, shown previously a pattern of disregulation caused by the treatment. The others show a pattern of KO overreaction to ROCKi
7.3.1 KO overreaction
Similar as with KO mimicking and KO rescue for identifying genes with KO pattern overreaction to ROCKi I used the three comparisons where KO treated is involved. (The comparison WT ROCKi vs KO ROCKi has the double weight as the others comparisons). Following this analysis I identify 225 genes with the pattern of overreaction
Interestingly only 19 genes are shred with the pattern of KO rescue and overreaction:
And this is the entichment of GO.
7.3.2 KO Unique-Response
With this comparison I want to highlight genes that reacted to ROCKi in WT but not in KO.
8 Rescue
\[ Rescue Ratio = \frac{log2FC(ROCKi_{KO})}{-log2FC(KO_{Veh})} \]
Where ROCKi_KO stands for the comparison of ROCKi vs Veh in KO genotype and KO_Veh is KO vs WT treated with Veh. A value around 1 (dashed line) represent complete rescue, values greater than 1 mean that not only WT Veh levels are restore but increased, values between 0 and 1 mean partial rescue, 0 means ROCKi did not altred KO expression, lastly negative values mean that with ROCKi the altered KO expression is increased with the same direction.
I talked with Nico about this and he pointed me to his paper where he and Ivano developed a Rescue Index which has this formula:
\[ RescueIndex=10^{-\frac{(Rescue - WT)}{(KO - WT)}} \]
Where normalized log2 counts were used. I decided to apply it to our data, but instead of power 10 I will use 2, for coherence with the powers. This is the formula:
\[ RescueIndex=2^{-\frac{(ROCKi_{KO} - Veh_{WT})}{(Veh_{KO} - Veh_{WT})}} \]
It is important to recall that log2(x)-log2(y)=log2(x/y), therefore I will use the log2-foldchanges calculated by DESeq2.
Here a value of 0.5 represent that there is no rescue and that the values from KO Veh and KO ROCKi change similarly in comparison to WT Veh. Values from 0-0.5 represent that the gene expression of KO ROCKi changed more than KO Veh (in the same direction). 1 means perfect rescue, KO ROCKi is similar to WT Veh, and values above 1 means over-rescue.
This is the relationship between the two values:
As it is clearly seen there is a relationship, where the data is only differentially transformed (maths not included, but analyzed). Here are the different groups plotted, and they were distributed as expected.
While writing this part just realized that I was not considering the actual effect of ROCKi in the formula so I came with a new one:
\[ Real\ Rescue\ Ratio = \frac{log2FC(ROCKi_{KO})-log2FC(ROCKi_{WT})}{-log2FC(KO_{Veh})} \]
The new term represent the difference between ROCKi and Veh in WT, the interpretation of this ratio is the same as the one I described above.
And the relation with Nico’s Rescue Index looks like this:
Patterns are maintained because as already discussed there are not that many and neither big expression changes between ROCKi and Veh in WT.
It is worth to mention that based on the density plot it is expected that most of the other genes show values of real rescue ratio similar as rescue genes:
Nevertheless, rescue pattern and no changes in gene expression would lead to similar Real Rescue Ratio. Here is an example, where both genes have a real rescue rate ~ 1.
Then, in order to make sense out of the Real Rescue Ratio I filtered for genes which expression needs to be rescued (in either direction) based on the log2 foldchange between KO and WT in Veh (abs(log2FC)>0.25).
With the real rescue ratio we could then visualize the results from Section 7.2.1. These are the genes with significant rescue pattern from the top 10 significant group-terms (sorted by p-value), the term with most genes with KO rescue pattern per group was selected as representative.
We can the plot all genes from a specific GO term and since there is a particular inetrest in the term “actomyosin structure organization” here it is (*yellow color represent genes that were excluded because they do not showed abs(log2FC) > 0. 25 in KO vs WT Veh):
We can additionally check KEGG pathways. Here are 3 examples, one is the Hippo signaling pathway and the others two where ROCK is involved:
Caution
Rho is colored yellow because KEGG assumes in this pathways only Rhoa is acting and the expression of Rhoa does not need to be rescued.
9 Model
Based on the results shown here I could infer this model
This might be too speculative.
10 KO animals
Talking with Fede about this she suggested to take a look at the mouse hippocampus data. As you know I could not find significant differential expressed genes between KO and WT mice.
10.1 GSEA
10.1.1 Stats
The first approach was to perform a GSEA on the “stats” from which p-value is calculated. Positives values would represent lower p-values. Extreme values represent lower p-values and the sign of of it, if gene is more expressed in KO (positive) or WT (negative). I tested here 6 lists, rescue genes, ROCKi insensitive KO genes altered and the merge of these two, for up and down patterns.
10.1.2 Fold Changes
Then I started testing independently up/down lists but here KO mice data used were log2 fold change between KO and WT. Positive values represent higher in KO than in WT.
WE HAVE SIGNIFICANT RESULTS!!!! The problem comes when comparing the sign of the enrichment score, If everything went perfect we would assume that downregulated genes would have been enriched in the downregulated genes of KO mice, and viceversa, but the opposite is happening.
The direct comparison is with the genes found differentially expressed ROCKi independant.
Stats correspond to pearson correlation.
R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=es_ES.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
time zone: Arctic/Longyearbyen
tzcode source: system (glibc)
attached base packages:
[1] grid stats4 stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] ggridges_0.5.4 ggraph_2.1.0
[3] igraph_1.5.0.1 rrvgo_1.12.0
[5] data.tree_1.0.0 IHW_1.28.0
[7] ggpubr_0.6.0 broom_1.0.5
[9] impute_1.74.1 ReactomePA_1.44.0
[11] ggsignif_0.6.4 ggforce_0.4.1
[13] circlize_0.4.15 gage_2.50.0
[15] GSVA_1.48.2 GO.db_3.17.0
[17] kableExtra_1.3.4 UpSetR_1.4.0
[19] enrichplot_1.20.0 pals_1.7
[21] patchwork_1.1.2 DT_0.28
[23] fgsea_1.26.0 pathview_1.40.0
[25] clusterProfiler_4.8.2 org.Mm.eg.db_3.17.0
[27] AnnotationDbi_1.62.2 data.table_1.14.8
[29] viridis_0.6.4 viridisLite_0.4.2
[31] ComplexHeatmap_2.16.0 cowplot_1.1.1
[33] EnhancedVolcano_1.18.0 DESeq2_1.40.2
[35] SummarizedExperiment_1.30.2 Biobase_2.60.0
[37] MatrixGenerics_1.12.2 matrixStats_1.0.0
[39] GenomicRanges_1.52.0 GenomeInfoDb_1.36.1
[41] IRanges_2.34.1 S4Vectors_0.38.1
[43] BiocGenerics_0.46.0 PCAtools_2.12.0
[45] ggrepel_0.9.3 lubridate_1.9.2
[47] forcats_1.0.0 stringr_1.5.0
[49] dplyr_1.1.2 purrr_1.0.1
[51] readr_2.1.4 tidyr_1.3.0
[53] tibble_3.2.1 ggplot2_3.4.2
[55] tidyverse_2.0.0 BiocManager_1.30.21.1
loaded via a namespace (and not attached):
[1] bitops_1.0-7 HDO.db_0.99.1
[3] httr_1.4.6 webshot_0.5.5
[5] RColorBrewer_1.1-3 doParallel_1.0.17
[7] numDeriv_2016.8-1.1 Rgraphviz_2.44.0
[9] tools_4.3.1 backports_1.4.1
[11] utf8_1.2.3 R6_2.5.1
[13] HDF5Array_1.28.1 mgcv_1.9-0
[15] lazyeval_0.2.2 apeglm_1.22.1
[17] rhdf5filters_1.12.1 GetoptLong_1.0.5
[19] withr_2.5.0 graphite_1.46.0
[21] gridExtra_2.3 fdrtool_1.2.17
[23] cli_3.6.1 scatterpie_0.2.1
[25] sass_0.4.7 labeling_0.4.2
[27] slam_0.1-50 KEGGgraph_1.60.0
[29] mvtnorm_1.2-2 tm_0.7-11
[31] askpass_1.1 systemfonts_1.0.4
[33] yulab.utils_0.0.6 gson_0.1.0
[35] DOSE_3.26.1 svglite_2.1.1
[37] dichromat_2.0-0.1 bbmle_1.0.25
[39] maps_3.4.1 rstudioapi_0.15.0
[41] RSQLite_2.3.1 treemap_2.4-4
[43] generics_0.1.3 gridGraphics_0.5-1
[45] shape_1.4.6 crosstalk_1.2.0
[47] zip_2.3.0 car_3.1-2
[49] Matrix_1.6-0 fansi_1.0.4
[51] abind_1.4-5 lifecycle_1.0.3
[53] yaml_2.3.7 carData_3.0-5
[55] rhdf5_2.44.0 qvalue_2.32.0
[57] blob_1.2.4 promises_1.2.0.1
[59] dqrng_0.3.0 bdsmatrix_1.3-6
[61] crayon_1.5.2 lattice_0.21-8
[63] beachmat_2.16.0 annotate_1.78.0
[65] KEGGREST_1.40.0 mapproj_1.2.11
[67] pillar_1.9.0 knitr_1.43
[69] rjson_0.2.21 codetools_0.2-19
[71] fastmatch_1.1-3 glue_1.6.2
[73] downloader_0.4 ggfun_0.1.1
[75] vctrs_0.6.3 png_0.1-8
[77] treeio_1.24.2 gtable_0.3.3
[79] emdbook_1.3.13 cachem_1.0.8
[81] openxlsx_4.2.5.2 xfun_0.39
[83] mime_0.12 S4Arrays_1.0.4
[85] tidygraph_1.2.3 coda_0.19-4
[87] pheatmap_1.0.12 SingleCellExperiment_1.22.0
[89] iterators_1.0.14 ellipsis_0.3.2
[91] nlme_3.1-162 ggtree_3.8.0
[93] bit64_4.0.5 bslib_0.5.0
[95] irlba_2.3.5.1 colorspace_2.1-0
[97] DBI_1.1.3 tidyselect_1.2.0
[99] bit_4.0.5 compiler_4.3.1
[101] rvest_1.0.3 graph_1.78.0
[103] NLP_0.2-1 xml2_1.3.5
[105] DelayedArray_0.26.6 shadowtext_0.1.2
[107] scales_1.2.1 rappdirs_0.3.3
[109] digest_0.6.33 rmarkdown_2.23
[111] XVector_0.40.0 htmltools_0.5.5
[113] pkgconfig_2.0.3 umap_0.2.10.0
[115] lpsymphony_1.28.1 sparseMatrixStats_1.12.2
[117] fastmap_1.1.1 rlang_1.1.1
[119] GlobalOptions_0.1.2 htmlwidgets_1.6.2
[121] shiny_1.7.4.1 DelayedMatrixStats_1.22.1
[123] jquerylib_0.1.4 farver_2.1.1
[125] jsonlite_1.8.7 BiocParallel_1.34.2
[127] GOSemSim_2.26.1 BiocSingular_1.16.0
[129] RCurl_1.98-1.12 magrittr_2.0.3
[131] GenomeInfoDbData_1.2.10 ggplotify_0.1.1
[133] wordcloud_2.6 Rhdf5lib_1.22.0
[135] munsell_0.5.0 Rcpp_1.0.11
[137] reticulate_1.30 ape_5.7-1
[139] stringi_1.7.12 zlibbioc_1.46.0
[141] MASS_7.3-60 plyr_1.8.8
[143] org.Hs.eg.db_3.17.0 parallel_4.3.1
[145] Biostrings_2.68.1 graphlayouts_1.0.0
[147] splines_4.3.1 hms_1.1.3
[149] locfit_1.5-9.8 reshape2_1.4.4
[151] ScaledMatrix_1.8.1 XML_3.99-0.14
[153] evaluate_0.21 httpuv_1.6.11
[155] tzdb_0.4.0 foreach_1.5.2
[157] tweenr_2.0.2 openssl_2.1.0
[159] polyclip_1.10-4 clue_0.3-64
[161] gridBase_0.4-7 rsvd_1.0.5
[163] xtable_1.8-4 reactome.db_1.84.0
[165] RSpectra_0.16-1 tidytree_0.4.4
[167] later_1.3.1 rstatix_0.7.2
[169] aplot_0.1.10 memoise_2.0.1
[171] cluster_2.1.4 timechange_0.2.0
[173] GSEABase_1.62.0